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Discrete dynamic models are a powerful tool for the understanding and modeling of large biological 
networks. Although a lot of progress has been made in developing analysis tools for these models, 
there is still a need to find approaches that can directly relate the network structure to its dynamics. 
Of special interest is identifying the stable patterns of activity, i.e., the attractors of the system. 
This is a problem for large networks, because the state space of the system increases exponentially 
with network size. In this work we present a novel network reduction approach that is based on 
finding network motifs that stabilize in a fixed state. Notably, we use a topological criterion to 
identify these motifs. Specifically, we find certain types of strongly connected components in a 
suitably expanded representation of the network. To test our method we apply it to a dynamic 
network model for a type of cytotoxic T cell cancer and to an ensemble of random Boolean networks 
of size up to 200. Our results show that our method goes beyond reducing the network and in most 
cases can actually predict the dynamical repertoire of the nodes (fixed states or oscillations) in the 
attractors of the system. 
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I. INTRODUCTION 

The interactions among cellular elements such as pro- 
teins, mRNA, and small molecules are orchestrated in 
such a way that they support all the complicated behav- 
iors cells are capable of (such as homeostasis, growth, 
movement, cell differentiation and cell division) p]. In 
order to get a full understanding of the relation between 
cellular behaviors and their underlying network of inter- 
actions, the construction of informative dynamic models 
based on the current biological knowledge is very impor- 
tant. Several dynamical modeling techniques exist, which 
provide different levels of detail in the dynamics, while in 
turn requiring varying amounts of biological information 
[213]. At one end of the spectrum, for example, highly 
quantitative information can be obtained from ordinary 
differential equation models [3HZ] by providing different 
reaction rates (e.g. transcription/translation rates, asso- 
ciation/dissociation constants, degradation coefficients, 
etc) and the biophysical/biochemical properties of the 
components (e.g. enzyme cooperativity) . At the other 
end, the qualitative dynamics of the system can be repro- 
duced by a discrete dynamic model ,8-13 , which requires 
only the combinatorial activating or inhibiting nature of 
the interactions, and not the kinetic details [Hj. 

Given the surprising but demonstrated fact that the es- 
sential dynamical properties of a variety of systems can 
be reproduced without knowing the values of the specific 
kinetic parameters of the processes involved [5l413j. one 
may wonder if there is a model-independent way to infer 
the dynamical properties of cellular networks just by us- 
ing the network topology (graph structure), that is, the 
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identity of the components and knowledge about their 
interactions. Historically, this relation between structure 
and dynamics was recognized early on in the pioneering 
work of Jacob and Monod [IS], Thomas (TB], Kauffman 
|17) . and Glass [IB], and was part of the original mo- 
tivation for the study of discrete dynamic models. The 
common idea is that the presence of feedback loops is nec- 
essary for the emergence of complex dynamical properties 
such as multistability and oscillations. More specifically, 
by assigning a sign to the interactions (-1- if activating and 
- if inhibitory) and to the feedback loops in the network 
(the sign of a loop is given by the product of the signs of 
its edges), the following two simple rules were proposed 
by R. Thomas ^9j to relate the network structure to its 
dynamics: 

1. A necessary condition for multistability (multiple 
stable steady states) is the existence of a positive 
feedback loop. 

2. A necessary condition for sustained oscillations 
(limit cycles) is the existence of a negative feedback 
loop. 

Since these early works, extensive research has been done 
in this direction and the validity of these rules has been 
demonstrated both in the differential [20H13] and dis- 
crete frameworks [2^26) . Recent works have even ex- 
tended these rules to include not only necessary but 
also sufficient conditions for multistability and oscilla- 
tions [571 dH]- 

Despite all this progress, there is still a need for devel- 
oping new analysis tools that relate the network struc- 
ture to its dynamics, especially ones that are applicable 
to large scale networks. This is a problem for many of the 
methods developed so far, since many of them are very 
computationally demanding and can only be exactly ap- 
plied to networks of small to moderate size. The size 



of the networks is also a problem even in cases where 
mathematical theorems are available, because as the net- 
work increases in size, it is very likely that the conditions 
needed in the theorems become harder and harder to be 
fulfilled. These limitations call for methods that are as 
generally applicable as possible. 

The novel analysis method we present in this work has 
the objective of inferring the dynamical repertoire of a 
network based purely on network topology and the com- 
binatorial nature of the interactions. Given a discrete 
dynamic model of a cellular network, the question that 
we try to answer is: what are the stable patterns of ac- 
tivity of the network (both fixed points and oscillations) ? 
Framed in the discrete dynamic framework, our method 
is based on the idea that some groups of nodes in the 
network can only stabilize in a small number of fixed 
states. By expanding the network to explicitly include 
the nature of the interactions (positive or negative) and 
the potentially synergistic regulation of every element in 
the network, we can identify these stable groups of nodes 
and use them to simplify the network. The result is a 
complete reduction (which directly gives the fixed points 
of the system) or a very simplified network in which most 
nodes are expected to oscillate. In the next sections we 
will explain our method in more detail, including the 
network expansion and network reduction techniques in- 
volved in it. 



II. PREDICTING THE STABLE DYNAMICAL 

REPERTOIRE OF A BOOLEAN MODEL OF A 

BIOLOGICAL NETWORK 



Biological networks and discrete dynamic 
models 



A network of cellular components can be repre- 
sented by a directed graph G = {V,E), where V = 
(f 1, U2, . . . , wjv) are the nodes describing the elements of 
the system, and E are the edges denoting the directed 
interactions among the components. To each edge one 
also commonly associates a sign, which denotes the reg- 
ulatory nature of the interaction (-1- if activating, and - if 
inhibitory). Although the sign of the interaction enriches 
the biological information included in the network, expe- 
rience shows that knowing the nature of the interactions 
is not enough and the combined effect of the interactions 
on each element also needs to be considered. For this 
purpose, every node Vi is assigned a function fi which de- 
pends on the ki regulators of Vi (and sometimes on itself) 
and that incorporates the combinatorial nature of the in- 
teractions. The set of functions F — {fi, f2, ■ ■ ■ , fw) con- 
tains all the dynamical information of the system, and 
thus, depends on the type of dynamic model used. 

For our study, we choose the simplest kind of dis- 
crete dynamic model, namely, the Boolean framework, 
in which the functions F are taken as logical (Boolean) 
functions and each node Vi can take one of two possible 



states: ON (or 1) and OFF (or 0). The biological inter- 
pretation of each state varies depending on the context, 
although in most cases it refers to above (ON) or below 
(OFF) a certain threshold level. The state of the system 
at any time can then be denoted by a vector whose i*'' 
component is the state of node Vi at that time. As a 
consequence of the discrete number of states, the state 
space of the system is finite (2^ states), and a temporal 
sequence of states can be represented as a trajectory in 
state space. Every temporal trajectory will eventually 
reach a set of network states in which it settles down, 
known as an attractor. An attractor can either be com- 
posed of a single network state in which the system stays 
fixed, known as a fixed point or a steady state, or a group 
of network states between which it alternates, usually re- 
ferred to as a complex attractor or an oscillation. In a 
steady state the state of all nodes remains fixed, while 
in a complex attractor all or a subset of the nodes keep 
changing their states (i.e. they oscillate), and the state 
of the rest of the nodes (if any) is fixed. 

The possible trajectories and attractors of a system 
depend not only on the functions F but also on the rep- 
resentation of time as a continuous or discrete variable. 
The most common choice for Boolean dynamics is taking 
time as a discrete variable, in which case the nodes are 
updated at discrete time steps according to the functions 
F. In the simplest case, the synchronous scheme, ev- 
ery node is updated simultaneously in discrete time steps 
and its state depends only on the state of the system in 
the previous step. The synchronous updating scheme, 
although suitable in some situations, is not apt for our 
purposes, as it inherently assumes that all processes oc- 
cur at a similar timescale, which is clearly not true for 
intracellular networks, in which a large variety of cellu- 
lar components and processes are involved. Furthermore, 
the timescales of a large number of these processes are not 
well known, and even when they are, they may be subject 
to fluctuations due to cell-to-cell variability and environ- 
mental perturbations. Previous work has developed var- 
ious asynchronous updating methods, wherein each node 
is updated according to its own timescale, and which can 
be either deterministic (the timescales are fixed during a 
simulation) j^H [S^ or stochastic (the timescales are ran- 
domly varied during a simulation) ^2S1 [5TI - I51] . In a com- 
parative study [35] several discrete time asynchronous 
updating schemes were tested in the same biological net- 
work, with its results suggesting that the general asyn- 
chronous method, in which at every discrete time step a 
randomly selected node is updated, is the most appropri- 
ate scheme for our requirements. 

The general asynchronous method is advantageous not 
only because it takes into consideration the multiple 
timescales involved in intracellular processes, our incom- 
plete knowledge of all these timescales, and the inher- 
ent stochasticity of biological processes, but also because 
it allows a natural biological interpretation of the at- 
tractors of the system. The reason for this is that, by 
definition, the general asynchronous method samples all 



possible timescales of the system, therefore, the attrac- 
tors must correspond to the patterns of activity of the 
system which are stable with respect to arbitrary fluctu- 
ations in the rates of the processes involved, which we 
identify as the stable dynamic repertoire of the network. 
The use of the discrete dynamic framework and the gen- 
eral asynchronous updating scheme then maps the prob- 
lem of finding the set of stable dynamic behaviors of a 
cellular network into finding the attractors of a discrete 
dynamic system. 



B. 



Finding the attractors of a discrete dynamic 
model 



As it was pointed out in section II A the state space of 
a Boolean network with N nodes contains 2^ states. This 
exponential dependence on the number of nodes of the 
state space's size makes the problem of finding the attrac- 
tors of a Boolean network computationally intractable, 
which means a full search of the state space can be per- 
formed, in practice, only for small networks {N < 20). 
To overcome this challenge several types of methods have 
been proposed to simplify the search space. The most 
prominent of these approaches, the so-called network re- 
duction methods [55H55] , shrink the effective network size 
by removing frozen nodes (nodes that reach the same 
steady state regardless of initial conditions) [38ti40] and 
dynamically irrelevant nodes (such as simple mediator 
nodes or nodes with no outputs), and by simplifying the 
Boolean functions (for example, due to the presence of a 
node with a fixed state), thus reducing the effective state 
space. 

Although the reduction methods developed so far have 
been successfully applied to several biological networks 
[551 - 1571 HT|. they have, however, some limitations. For 
example, some of them may not be effective enough, in 
the sense that even after reduction the state space's size 
is still unmanageable, in which case one must resort to 
sampling the state space [551 SI]- Other methods can 
affect the state space in such a way that the attractor 
space is changed and only some types of attractors (e.g. 
steady states) are preserved [5S1 [57] . 



C. The role of stable motifs and oscillating 
components in the attractor landscape 

The method that we propose is based on the idea that 
certain components of the network can only stabilize in 
one or a small number of attractors. This idea is itself 
not new and is closely related to R. Thomas' rules linking 
feedback loops with the appearance of complex dynam- 
ical behavior in biological networks [H]. For example, 
an approach that connects the dynamics of certain net- 
work motifs to construct the attractors of the full Boolean 
network was recently proposed by Siebert |28| . The nov- 
elty of our method is the efficiency of identifying every 



motif that stabilizes in an attractor despite being cou- 
pled to the rest of the network. The main insight of our 
approach is that a representation of a Boolean network 
known as the expanded network can be used to easily 
identify these network components and the states they 
can stabilize on in an attractor. By combining the knowl- 
edge of the stable behavior of this group of nodes with 
network reduction methods, we can find other network 
components that stabilize as a consequence of the for- 
mer. By repeatedly applying this procedure, the stable 
states of all nodes can be found, which will correspond 
to their states in the attractors of the system. 

For ease of understanding we will focus on finding the 
components that stabilize in a fixed state. More specif- 
ically, we will look for network components with cer- 
tain topological characteristics that cause themselves and 
other nodes to take a fixed state. We refer to these net- 
work components as stable motifs or stable components. 
It is worth noting that the nodes in the stable motifs 
may or may not be part of the so-called frozen nodes of 
Boolean networks [551110] , since the nodes of these stable 
motifs can have more than one steady state. Although 
it may seem at first that this restriction to fixed state 
nodes reduces the general applicability of the method, 
it turns out that this is not necessarily the case. In- 
deed, there are two possibilities after using our method 
to find all the nodes that take a fixed state in an at- 
tractor. The list of fixed-state nodes could include all 
nodes in the network, in which case we have identified 
a fixed point attractor. Or, if the list covers a subset of 
the nodes, these nodes must represent the non-oscillating 
nodes of a complex attractor. In this case the nodes the 
method doesn't identify are expected to keep changing 
their values (i.e. oscillate) in the attractor. Analogously 
to the stable case, we refer to the network components 
that stabilize in an oscillating state as oscillating motifs 
or oscillating components. We discuss in more detail the 
role of these oscillating components and of oscillations in 
section IIIGI 



D. Expanded Network 

In order to identify the stable motifs of a Boolean net- 
work, it is convenient to use a representation that incor- 
porates explicitly the update functions F. Previous work 
[^j has shown that a useful representation for this pur- 
pose is the so-called expanded network representation. 

The creation of the expanded network consists of two 
basic operations, which we illustrate in Figure[l] First, to 
explicitly take into account negative regulations, we in- 
troduce a complementary node Vi for every node Vi in the 
network and assign to each Vi an update function which 
is the Boolean negation of Vi 's update function fi . The 
addition of complementary nodes has a two-fold effect; 
not only does it allow us to evaluate the inhibitory effect 
of a node on the rest of the network, but by assigning the 
negation of the original update function to every comple- 



(a) 

!a = NOT B 



(b) 



(c) 



♦ 



fA = B 
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FIG. 1: Operations for the creation of the expanded network, (a) Node A has the update function /a = NOT B. 

The addition of complementary nodes introduces a new node A with update function /^ ~ NOT Ja — B. It also 

introduces a complementary node for node i3, which makes the update function of A take the form 

Ja = NOT B = B. (b) Node A has the update function fA = B AND C. The addition of composite nodes 

introduces a new node BC that represents the cooperative effect of B and C on A. (c) Node A has the update 

function Ja = B AND (NOT C). The two expansion operations introduce complementary nodes for A, B, and C, 

and a composite node for the AND relation between B and C. 



mentary node, it also explicitly considers how the other 
nodes can promote the inactivation of a given node. Note 
that because of the addition of complementary nodes all 
the edges are of the same (positive) nature, and thus, no 
sign needs to be specified. 

As an example, let us consider node A with update 
function Ja = NOT B, as illustrated in Figure Ilia), 
where for simplicity we denote the state of the node with 
the node name. The addition of complementary nodes 
means that we add a new node A with update function 
/j = NOT fA = B. An additional complementary node 
is added for node B, which makes the update function 
of A take the form /a = NOT B = B. Th_e expanded 
network contains two positive edges, from B to A, and 
from B to A, instead of the negative edge from B to A. 

Second, to incorporate the combinatorial nature of the 
update functions, we introduce a composite node for each 
set of synergistic interactions (that is, AND relation- 
ships) in the Boolean functions F. For example, consider 
the case shown in Figure [lib) , in which node A has the 
logical function fA = B AND C. Since the function con- 
tains an AND relationship between node B and node C, 
a composite node BC is added when expanding the net- 



work. Node B and C are connected by directed edges to 
the composite nodes BC, and an edge from BC to A is 
also added. A more complicated example in which both 
operations are applied is shown in Figure [lie) . 

In general the introduction of composite nodes may 
not be as obvious as shown in these examples if we have 
nontrivial combination of AND, NOT and OR rules in 
the functions F. Because of this, it is convenient to 
represent each update function fi with K input nodes 
{viiiVi2i ■ ■ ■ -iVix} in the following disjunctive normal 
form: 

/ = (si AND S2 AND • • • AND sj) 
OR (sj+i AND Sj+2 AND ■ ■ ■ Sk) 
OR • • • OR (s; AND si+i AND • • • AND s„) , 

where the s^'s are either the states of one of the K in- 
put nodes of /, or one of these states' negations. In the 
same way, the negation of the update functions, /j, is 
also represented in a disjunctive normal form. Once the 
functions are represented in the disjunctive normal form, 
the introduction of composite nodes is simple: a com- 
posite node is added for every set of nodes involved in 



an AND-dependent relationship. In the foUowing we wiU 
refer to nodes that are not complementary nor composite 
as normal nodes. 



E. Identifying Stable Motifs from the Expanded 

Network 



We define a stable motif in the expanded network 
as any of the smallest strongly connected components 
(SCCs) in the expanded network representation which 
satisfy these two properties: (1) the SCC does not con- 
tain both a node and its complementary node, and (2) 
if the SCC contains a composite node, all of its input 
nodes must also be part of the SCC. The first condition 
makes sure that there is no contradiction between the 
SCCs found and a state in the original Boolean network, 
wherein every node can either take the value (which 
would correspond to having the complementary node in 
the SCC) or 1 (which would correspond to having the 
normal node in the SCC). The second condition is a con- 
sequence of the synergistic nature of composite nodes, 
which means that a composite node and all of its inputs 
form an irreducible unit. By smallest SCC we refer to 
any SCC that does not contain another SCC with the two 
specified properties, but that, otherwise, is arbitrary in 
size. For example, a node of the expanded network with 
a self-loop would be one of these smallest SCCs. We re- 
strict ourselves to the smallest SCCs so that when there 
are multiple such SCCs in the system all possible combi- 
nations of these SCCs are considered as steady states of 
the network. Note that this docs not result in the loss of 
information on larger possible SCCs; the remaining parts 
of these SCCs will be found after the smallest SCCs are 
reduced. 

The composition of the stable motif directly deter- 
mines the state of a subset of nodes in the Boolean net- 
work: every normal node of the stable motif will adopt 
the state 1, and for every complementary nodes included 
in the stable motif the corresponding Boolean node will 
adopt the state 0. 

As an example, consider the Boolean network and its 
expanded network representation in Figure [21 The small- 
est SCCs that satisfy both stable motif requirements are 
JA, S, C} and {B, C}, both of which are shown in Figure 
[2Fc) . The corresponding states for these stable motifs are 
JA^1,B = 1,C = 0} a.nd{B = 0,C = 1}, respectively. 

So far, we have only defined a stable motif in terms 
of the expanded network representation. We can extend 
the concept of stable motif to the original Boolean net- 
work to mean the nodes in the Boolean network whose 
state is specified by a stable motif of the expanded net- 
work. These nodes include all the normal nodes that are 
included in the stable motif of the expanded network and 
all the normal nodes whose complementary nodes are in- 
cluded in the stable motif of the expanded network. In 
this way, a stable motif can mean either a set of nodes in 
the expanded network that satisfy the two requirements 



or the set of nodes in the Boolean network whose state 
is specified by a stable motif in the expanded network, 
depending on the context. 

It is important to point out that the stable motifs de- 
pend on the structure of the logical rules and thus on the 
topology of the network being considered. This means 
that an arbitrary change in the logical rules or in the 
topology of the network can modify the stable motifs, 
and we know of no obvious way to determine how the 
motifs will change without having to reconstruct the ex- 
panded network. 



F. Netw^ork Reduction 

Once the stable motifs of the network have been iden- 
tified, the next step is to determine the infiuence of these 
nodes on the rest of the network. More specifically, for 
each stable motif found, we want to find the nodes in 
the network whose state stabilizes due to the influence of 
this stable component. We adapt the method previously 
developed by Saadatpour et al. to simplify the network 
[35l 141] , which has been shown to preserve both the fixed 
points [351 137] and the complex attractors of the system 
|56) . This method removes not only the frozen nodes of 
network [38ti40] . but also the nodes that reach a steady 
state under the influence of a given combination of source 
node states. It consists of two main steps: 

1. Identify the nodes whose value is fixed during the 
dynamics, which we will refer to as source nodes; for 
our case these will initially correspond to the nodes 
in the stable motif being considered. Modify the 
Boolean functions of the nodes downstream of the 
source nodes by setting the value of the source node 
to its fixed value. If any of these modified node's 
functions can only have one possible outcome, then 
this node can be used as a source node itself. 

2. Remove mediator nodes (i.e., nodes that have only 
one incoming edge and one outgoing edge) and ir- 
relevant sink nodes (i.e., nodes that have no outgo- 
ing edges) . For the case of mediator nodes connect 
the input of the mediator node to its output. The 
value of these removed nodes will be determined 
once the value of their input nodes is known. 

For each separate stable motif found in the expanded 
network, these two steps are repeated recursively until 
neither of them can be applied anymore. 

After network reduction, we obtain a set of states for 
each stable component, each of which corresponds to the 
states of the nodes in the stable motif and the states of 
other nodes which stabilized as a consequence of the sta- 
ble motif. For each of these sets of states, there is also 
a reduced network that contains the nodes whose value 
we still do not know. On each of these reduced networks 
the whole method will be applied again, starting with the 
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FIG. 2: Identification of stable motifs from the expanded network, (a) An example of a Boolean network, (b) The 

expanded network representation of the Boolean network in (a), (c) The two stable motifs in the expanded network, 

that is, the two smallest SCCs in the network that satisfy the requirements of not containing both a node and its 

complementary node, and containing all the inputs of every included composite node. Each stable motif indicates 

the fixed states of the corresponding subset of Boolean nodes. 



creation of the expanded network (section II D I and end 



ing with the network reduction process, and iteratively 
until there are no more nodes with unknown states or no 
new stable motifs are found. For the case where there 
are no more nodes with unknown states, a fixed point at- 
tractor of the system is obtained directly from the state 
of the nodes of the stable components. 

For the cases in which there are no new stable mo- 
tifs in the final reduced networks, the state of the nodes 
making up said networks are still unknown. Since our 
method is based on identifying nodes that stabilize in a 
specific steady state, the expectation is that these left- 
over nodes will oscillate in an attractor of the system, 
while in that same attractor the rest of the nodes will 
take the steady state value found during the simplifica- 
tion process that leads to the reduced network in consid- 
eration. For conciseness we will refer to the final output 
of our method, consisting of a set of stabilized nodes (and 
their states) and a (potentially empty) set of nodes with 
undetermined states as a quasi- attractor. Our notion of 
quasi-attractor is closely related to similar other concepts 
in the literature, such as the singular steady state origi- 
nally introduced by Snoussi and Thomas [128, 43! and the 
logical steady state used by Klamt et al. |44j . 



Quasi-attractors are closely related to the attractors of 
network, both fixed points and oscillations. For example, 
if the set of oscillating nodes in a quasi-attractor is empty, 
then the states of the stabilized nodes will correspond 
to the node states in a fixed point attractor, thus, this 
quasi-attractor is in fact a fixed point. More generally, 
for every attractor of the system there exists a quasi- 
attractor associated to it; this quasi-attractor is such that 
every node whose state is fixed in the quasi-attractor will 
also have its state fixed in the same value in the attractor 
it is associated to. The proof for this is statement is 
shown in Appendix \K\ 



G. Oscillations and oscillating components 

Using the expanded network representation on net- 
works that show oscillatory behavior, we have found that 
it can also be used to identify the oscillating components 
of a network. To find the oscillating components using 
the expanded network representation, we search for the 
largest SCCs that satisfy these properties: (1) the SCC 
must contain the complementary node of every normal 
node and viceversa, and (2) if the SCC contains a com- 



posite node, all its input nodes must also be part of the 
sec. The first of these conditions makes sure that all 
nodes oscillate, by having both of the states of every node 
as part of the SCC. The second condition is a consequence 
of a composite node and all of its inputs forming an ir- 
reducible unit. In this case we look for the largest SCCs 
because we want to find all the nodes that feedback to 
each other in the oscillation. 

These properties are necessary but not sufficient con- 
ditions for a group of nodes to oscillate. We have found 
that there is a third condition that, if also satisfied, is 
sufficient (though not necessary) for a group of nodes to 
oscillate, which is that (3) there can be no stable mo- 
tifs that contains no composite nodes as part of the os- 
cillating component. This extra condition is related to 
the possibility of the coexistence of a steady state and a 
complex attractor in the sub-state-space. The simplest 
example that shows this kind of behavior, which we term 
unstable oscillations, is shown in Figure [3J In general, 
during the reduction process, we need to find the compo- 
nents that could have unstable oscillations (that is, that 
satisfy the sufficient conditions (1) and (2), but not the 
necessary condition (3)) to make sure that we preserve 
all attr actors. 

Another type of dynamical behavior of the oscillat- 
ing components that needs to be considered is when the 
nodes of the oscillating components do not visit all possi- 
ble states of their sub-state-space in an attractor, which 
we refer to as an incomplete oscillation. Incomplete oscil- 
lations are important because a node that is downstream 
of an oscillating component that displays incomplete os- 
cillations may reach a steady state as a consequence of 
the nodes of the component only visiting part of their 
sub-state-space. This type of behavior has been found 
before in studies of synchronous networks (for example, 
see Figure 1 in the work by Bilke and Sjunnesson [38 J. 

As an example of an incomplete oscillation, consider 
the network shown in Figure pta) . In this example the 
nodes A and B oscillate and their state transition graph is 
shown in Figure Efb). From their state transition graph 
one can clearly see their complex attractor is A, B — 
{(1, 0), (0, 0), (0, 1)}. Since once the nodes A and B settle 
down in the attractor, the state A = \, B = \ cannot be 
reached, then node C (whose update function is jc — 
A AND B) will necessarily stabilize in the state C — 0. 
Note that if either A and B took all possible states in the 
attractor, or if the update function of C was different, C 
would also oscillate in the attractor. 



III. RESULTS 

We implement the network reduction process with a 
custom Java code. The steps of the network reduction al- 
gorithm are described in AppendixlB] The main challenge 
in implementing the reduction method computationally 
lies in finding the stable motifs from the expanded net- 
work representation. The reason for this is that stable 



motifs have to be the smallest SCCs that satisfy the prop- 
erties we outlined above, which means that in order to 
identify all possible stable motifs we need to find all pos- 
sible directed cycles, since each of them could potentially 
be the smallest SCC we are searching for. The issue with 
finding all possible directed cycles is that the time com- 
plexity is O {{N + E){c + 1)) (using Johnson's algorithm 
|46|). where N is the number of nodes, E is the number 
of edges, and c is the number of directed cycles, the lat- 
ter of which can grow faster than 2^ for the worst case 
scenario of a fully connected network. 

Because of the caveats discussed in section HTG] involv- 
ing unstable and incomplete oscillations, one may be con- 
cerned that other similar cases are not taken into account 
and that, as a consequence, some attractors could be lost 
during the reduction process. In Appendix \K\ we address 
this concern by formally proving that for every attractor 
in the Boolean network there is a corresponding quasi- 
attractor that will be found by our reduction method. In 
order to further test the validity and generality of our 
network simplification method, we apply it to a previ- 
ously developed biological cellular network, and also to 
an ensemble of random Boolean networks [T71 |3S] . For 
the case of the biological network, we choose the signaling 
and regulatory network involved in a type of white blood 
cell cancer (T cell large granular lymphocyte leukemia or 
T-LGL leukemia) |TTJ |5T] , while for the ensemble of ran- 
dom Boolean networks we choose the original Kauffman 
ov N - K model fT7l|45l. 



T Cell Large Granular Lymphocyte Leukemia 
Network 



Cytotoxic T cell lymphocytes (CTLs) play a central 
role in the immune response. When an infection occurs, 
CTLs detect antigens in infected cells and, in response, 
trigger a set of intracellular signaling cascades, which lead 
to the production of cytokines (small signaling molecules) 
that induce the self-destruction of the infected cells. Nor- 
mal cytotoxic T cells undergo activation-induced cell 
death (or apoptosis) after successfully fighting infection, 
however, in T-cell large granular lymphocyte (T-LGL) 
leukemia mature cytotoxic T lymphocytes survive and, 
in time, cause an autoimmune disorder. In addition to 
their abnormal survival, these CTLs also show a deregu- 
lated activity (higher or lower than in normal CTLs) of 
many signaling pathways and genes. 

A Boolean network model of T cell survival in the con- 
text of T-LGL leukemia was constructed by Zhang et al. 
[TT] through an extensive literature search. The logical 
rules were constructed such that the known experimen- 
tal results in healthy and leukemic cytotoxic T cells were 
reproduced by the model and in many cases do not have 
a simple form (the logical rules are reproduced in Ap- 
pendixO . The resulting network consists of 60 nodes and 
142 regulatory edges, with the nodes representing genes, 
proteins, receptors, or small molecules (Figure pi. The 
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FIG. 3: An example of a component that has unstable oscillations. This network has an attractor in which all the 

nodes oscillate and also has a steady state attractor. (a) The network and its respective Boolean rules, (b) The 

state transition graph of the network. The nodes of the state transition graph are the states of the system (written 

in the order A,B) and the edges represent allowed state transitions when only one node is updated. State 11 is a 

fixed point as there are no transitions going out of it. States 01, 00, and 10 form a complex attractor. (c) The 
expanded representation of the network. Note that {A, B, AB} forms a stable motif and that the whole expanded 

network forms an oscillating SCC. 



(a) 




(b) 



C 



01 M- 



t) 



11 



i 




/a = (NOT A) OR (NOT B) 
Jb = (NOT A) OR (NOT B) 
fc = A AND B 

FIG. 4: An example of a node configuration in which a node can stabilize without the influence of an input signal or 

a stable motif. In this example A and B oscillate in a complex attractor, but they do not take all possible states of 

their state transition graph in this attractor. Specifically they miss the A = 1, B = 1 state. As a consequence node 

C stabilizes in the state C = 0. (a) The node configuration and their respective Boolean rules, (b) The state 

transition graph of nodes A and B. States 01, 00, and 10 form a complex attractor. 



network contains 6 nodes with no upstream components 
which represent external signals (Stimuli, IL15, PDGF, 
Stimuli2, CD45, and TAX), and also contains 3 output 
nodes that serve as indicators of biological functions or 
cell fate (Cytoskeleton signaling. Proliferation and Apop- 
tosis) . Two of these input and output nodes play an es- 
pecially important biological role: Stimuli, which repre- 
sents antigen stimulation, and Apoptosis, which denotes 
programmed cell death. 

Zhang et al. [TT] used asynchronous Boolean dynamics 
to show that in the sustained presence of PDGF and IL15 
the system may converge to a state that recapitulates all 
dysregulations in T-LGL leukemia, in addition to the ex- 
pected state of self-programmed cell death (apoptosis). 



Later Saadatpour et al. [55] used network reduction to 
show that, under the presence of said signals, the two 
attractors found by Zhang et al. (apoptosis and T-LGL 
leukemia) are the only possible ones. Although in the 
case studied by Saadatpour et al. the network reduction 
method was enough to simplify the network to a man- 
ageable size (6 nodes), this is actually not the case if 
one is interested in studying all possible combinations of 
the input signals, since in many cases the network ob- 
tained after reduction is still quite large (30-40 nodes). 
This then gives an opportunity to apply our reduction 
method to cases in which previous methods fall short. 

We apply our reduction method to all combinations 
of external signals in the presence of antigen (Stim- 
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FIG. 5: The T-LGL leukemia survival signaling network. The shape of the nodes indicates the cellular location or 
the type of nodes: rectangles indicate intracellular components, ellipses indicate extracellular components, diamonds 

indicate receptors, and hexagons represent conceptual nodes (Stimuli, Stimuli2, P2, Cytoskeleton signaling. 

Proliferation, and Apoptosis). Node colors are used to distinguish input nodes (white), output nodes (black) and the 

rest of the nodes in the network (gray) . An arrowhead or a short perpendicular bar at the end of an edge indicates 

activation or inhibition, respectively. This figure and its caption are adapted from |41| . 



uli=ON). To validate the quasi-attractors found through 
our method, we compare them with the attractors ob- 
tained by randomly sampling a large number of ini- 
tial conditions and evolving them using a general asyn- 
chronous updating scheme wherein one node is updates 
at each time step. For all the cases we find that 



the attractors/quasi-attractors obtained are exactly the 
same, which together with the proof in Appendix [^ 
strongly suggests that the reduction method can indeed 
be used to find all possible attractors. A table contain- 
ing all the leukemic attractors is included in Appendix [D| 
(the Apoptosis=ON attractor, in which all nodes except 



10 



Apoptosis are inactive, is always a possibility, thus for 
simplicity we do not include it in the table). 

As an example, consider the case in 
which the IL15 signal is present and the 
rest of them are not (IL15=Stimuli=0N, 
PDGF=Stimuli2=:CD45=TAX=OFF) . Simplifying 

the network using only the effect of these input sig- 
nals we obtain a Boolean network of 42 nodes. The 
expanded network representation of this network has 
144 nodes (42 normal nodes, 42 complementary nodes, 
and 60 composite nodes) and 302 edges. Searching 
the expanded network, we find four stable motifs from 
the expanded network representation which correspond 
to the states: i) {PDGFR = SIP = SPHKl = ON, 
Ceramide^OFF}, ii) {PDGFR = SIP = SPHKl = 
OFF} , iii) {TBET = ON}, and iv) {P2 = ON}. 
These motifs are shown in Fig. [61 Performing network 
reduction using each of these four stable motifs leads to 
reduced networks of widely varying sizes: 3, 39, 35, and 
39 nodes, respectively. The reduced network due to the 
first stable motif consists of three disconnected nodes 
with self-loops (one negative and two positive ones), and 
can gives rise both to apoptosis or the T-LGL leukemia 
attractor. For the networks corresponding to the three 
remaining motifs we need to continue the reduction 
process and search for stable motifs in each of these 
networks. 

The stable motifs we find during reduction of the re- 
maining networks are likely to include some of the same 
four stable motifs found previously (as long as that spe- 
cific motif was not used to obtain the specific network 
in consideration) , but may also contain stable motifs dif- 
ferent from those previously found. For example, the 
network due to the second motif has two stable motifs, 
both of which were found in the previous network (iii and 
iv) . The network obtained from the third motif has four 
stable motifs, one of which is different from the motifs 
previously found, with states ({MEK = ERK = RAS 
= PCLGl = IL2RBT = IL2RB = GRB2 = ON}), and 
three of which are the same as previous motifs (i,ii, and 
iv). For the network corresponding to the fourth motif 
we find three stable motifs, all of which had already been 
found (i, ii, and iii). If we continue the reduction process 
we find that the network due to the second motif gives 
rise to the apoptosis attractor after 2-3 more network 
reductions (depending on the stable motifs used for the 
reduction), while the third and fourth motifs can produce 
both the T-LGL leukemia or apoptosis attractor after 2-4 
more network reductions. 



B. Ensemble of Random Boolean Netw^orks 

Random Boolean networks were first introduced by S. 
Kauffman as a model to understand the general dynam- 
ical properties of gene regulation and cell differentiation 
processes [T7], and have been extensively studied ever 
since [JS]. In addition to the original Kauffman net- 



works, several variants of random Boolean networks that 
could be considered more biologically realistic have also 
been introduced (for example, models with arbitrary de- 
gree distributions [37], canalizing functions [3H], thresh- 
old functions \^\, or multiple discrete states [SD]). All 
of these models share the distinguishing feature of the 
original model, that is, the existence of three dynamical 
regimes: i) an ordered one in which similar initial condi- 
tions typically converge after a transient time, ii) a disor- 
dered regime in which the system becomes very sensitive 
to small changes in the initial conditions, and iii) a criti- 
cal regime, poised at the boundary of the ordered and dis- 
ordered regimes, in which perturbations retain their size. 
Evidence suggests that the gene regulatory networks of 
living organisms operate near the critical regime [511 I52j . 

For simplicity we use the original Kauffman or N — K 
model to test our network reduction method in randomly 
constructed networks. The N ~ K model consists of an 
ensemble of Boolean networks with N nodes in which ev- 
ery element has K input nodes. To construct one of the 
networks in this ensemble, the K input nodes of every 
element are chosen randomly from the rest of the net- 
work. Every node is then assigned one of the 2^ possible 
Boolean functions at random. To use what is considered 
the most biologically realistic case of this model, we use 
different network sizes with degree K = 2, which is the 
case at which this ensemble operates in the critical regime 
|45) . For this case it has been shown that the number of 
relevant nodes increases as N^^^, and that the number 
of asynchronous attractors grows as a power law in N 
|58) : one would then expect the existence of an efficient 
method to find this relatively small number of attractors 
with a relatively small fraction of relevant nodes. 

To test the validity of our network reduction method, 
we compare the final reduced networks obtained with the 
asynchronous attractors of the original network for an en- 
semble of Kauffman networks of different sizes {N ~ 5, 
10, 15, 18, 25, 50, 100, 150 and 200, with an ensemble 
size oifl = 200 for iV < 18 and fi = 100 for N ^ 25). 
For networks up to size N = 100 we were able to use 
the exact Johnson's algorithm to find the networks' cy- 
cles, however, for larger networks we needed to restrict 
the search to cycles of less than 40 nodes. To compare 
the system's attractors with the result of our reduction 
method (the quasi-attractors) we focus only on the nodes 
whose state stabilizes; if for every quasi-attractor there is 
one attractor that exactly matches the stabilized states 
of the quasi-attractor, we then say the results are com- 
patible. Note that since our reduction method does not 
predict the actual state of the nodes remaining after re- 
duction of the stable motifs, it is not necessary for these 
node states to agree. The naive expectation is that these 
remaining nodes oscillate in the attractors. If this is in- 
deed the case, we then say that the results are equivalent. 

To find the attractors for small networks (TV < 18), we 
construct the asynchronous state transition graph, a di- 
rected graph on the unit TV-cube whose nodes represent 
the states of the system and whose edgess are the allowed 
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FIG. 6: The three stable motifs of the T-LGL leukemia network found most often during the reduction process. The 

actual motifs found and the states in which each of these motifs can stabiUze vary depending on the active signals. 

We also show the input signals (white nodes) that affect these motifs directly or almost directly (for the motif in 

(c)). (a) The PDGFR-SlP-SPHKl-Ceramide motif, which represents the ceramide/sphingomyelin pathway and the 

platelet derived growth factor receptor, (b) The IFNG-P2 motif, which is related to the control of the cytokine 
interferon gamma in CTLs. (c) The TBET motif, which represents the regulation of the T-box transcription factor. 



transitions between states as a result of a single node's 
update [SI] [32]. In the asynchronous state transition 
graph attractors correspond to sink SCCs, that is, SCCs 
of states whose outgoing edges can only lead to other 
states of this same SCC. More specifically, fixed points 
correspond to single states with self-edges and no other 
outgoing edges, while complex attractors correspond to 
sink SCCs made up of more than one state. 

For large networks {N > 25) we cannot construct the 
whole asynchronous state transition graph because of its 
enormous size, so we resort to sampling the state space 
to look for attractors. In particular we use a method 
based on the one by Wang et al. [57]. In this method 
we construct part of the state space by starting from a 
large number Ng of initial conditions and following the 
system's trajectory for T effective time steps, that is, we 
make sure that at every step one node changes value (un- 
less a fixed point is reached, in which case no nodes can 
change value). To find the attractors from the resulting 
partial state transition graph we use the same criteria as 
in the complete state transition graph. To avoid false 
positives, we check the validity of every attractor ob- 
tained with this method by starting from one of the states 
in the putative attractor, updating it Ttransient effective 
time steps, then creating a partial state transition graph 
with Tsearch effective time steps, and finally searching for 
attractors in the resulting state transition graph. For our 
case we use Ns = 5000, T — 300, Ttransient = 1500 and 

Tsearch = 50000. 

Of the total 1300 networks that were compared we find 
that in all but five networks (all of which had iV > 100) 
the results of our method and of the attractor identifica- 
tion/sampling methods were equivalent. For the remain- 
ing five networks we find that the results were compati- 
ble, that is, although the state of the nodes predicted to 
have stabilized by the reduction process matched in the 
results of both methods, there were some nodes that did 
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not stabilize in the reduction method that were found to 
be stable in the attractors found by sampling (i.e. using 
the partial state transition graph) . We reiterate that this 
disagreement does not mean that our method is incorrect; 
our reduction method does not actually predict the state 
of the nodes remaining in the final reduced networks. It is 
also worth pointing out that in most quasi-attractors the 
fraction of nodes that don't stabilize is relatively small, 
as shown in Figure [Tj Although these remaining nodes 
are expected to oscillate in the attractor, this is actually 
not necessary, as discussed in section [IIG[ 

We also compare the number of attractors found by 
each of the three methods. For small networks we find 
that network reduction and the exact methods always 
find the same number of attractors/quasi-attractors. For 
large networks we find that the reduction method al- 
ways finds either more or the same number of attrac- 
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than the sampling method. For all network sizes shown 

an ensemble size of 100 networks was used. 



tors than the sampling method. In Figure |8] we show 
the average, the 20th percentile, and the 80th percentile 
of the difference in the number of attractors found by 
the two methods for large networks. To make sure that 
the quasi-attractors found by reduction are real attrac- 
tors, we check their validity by constructing a partial 
state transition graph just as we did for the sampling 
method. In some cases, in which the putative attractor 
was expected to have a large number of oscillating nodes 
(>20), an attractor could not always be found with this 
method. For these cases we analyzed the trajectory in the 
partial state transition graph and identified which nodes 
changed states and which of them didn't, and compared 
them with the putative attractors. In all cases we find 
that the results are equivalent (or compatible for the five 
networks mentioned before). This result, together with 
our proof in Appendix |Xj strongly suggests that each 
quasi-attractor of the reduction method correspond to 
an attractor of the system. 

Finally we compare the time performance of the meth- 
ods. In Figure ^a) we show the average time to comple- 
tion of the three methods on the same ensemble for differ- 
ent network sizes. Although at very small network sizes 
{N < 10) the exact method (the whole state transition 
graph) is on average faster than the reduction method, 
for larger networks the reduction method is faster than 
the others. For large networks the reduction method is, 
not only faster on average than the sampling method, 
but the distribution of times shows that, for all sizes, 
network reduction takes less than a second for most of 



the networks (see Figure ^b) for the N — 100 case). 
This is true even at larger network sizes. For example, 
for N = 100, 150, and 200 we have 88%, 76%, and 71% of 
the networks, respectively, take less than one second. In 
contrast, for the sampling method none of the networks 
of these same sizes take less than one second. These re- 
sults suggest that our method is not only more effective 
than state space sampling in the sense that it does not 
miss any of the attractors in the system, but it is also 
significantly more efficient in terms of time performance. 



IV. DISCUSSION 

In this work we have presented a novel reduction 
method that can greatly simplify a network, allowing us 
to deal with system sizes of an order of magnitude larger 
than what is possible through full state space searching 
methods. This reduction method uses an expanded rep- 
resentation of the network that explicitly includes the 
nature and logic of the interactions, which allows us to 
identify the network motifs that stabilize in a steady state 
and use them to simplify the network. To test the validity 
of our method we applied it to a biological network (the 
T-LGL leukemia survival network) and an ensemble of 
random Boolean networks of various sizes. We find that 
the results of our method always agree with the behavior 
of all networks tested. 

An important point to make is that our reduction 
method does not actually reduce the complexity of the 
attractor-finding problem. What it does is to transfer 
the complexity of the problem from the state space size 
to how complex the network is, specifically, to how many 
cycles the network contains. The reasoning behind this 
transfer is that we want to take advantage of the sparse- 
ness of biological networks to make the problem more 
tractable for larger network sizes than a brute force ap- 
proach (i.e. sampling the whole state space) would ac- 
complish. Indeed, our results show that our reduction 
method is able to outperform a sampling of the state 
space both in terms of attractors/quasi-attractors found 
and in terms of computational time (see Figures p^ and 
[9]). Added to this is the fact that in not a single case 
did state space sampling find attractors not compatible 
with the results of network reduction, as was expected 
from our proof that the reduction method preserves all 
attractors. 

A surprising result of applying our method to the en- 
semble of random networks was that for most cases it ac- 
tually predicted the full dynamic repertoire of these net- 
works, that is, it found which nodes stabilize and which 
of them oscillate in the attractors. This agreement is not 
trivial because our reduction method does not predict 
the state of the nodes remaining in the final reduced net- 
works. What is the reason for this success? To answer 
this question, we need to remember that our reduction 
method is based on finding stable motifs and using them 
to simplify the network. Therefore, in the final reduced 
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FIG. 9: Time performance of the different methods (see also the main text), (a) The average time it takes to find 
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N = 150 case for the reduction method is the consequence of a network in the ensemble that took an unusually long 

time because of the large number of cycles in the network, (b) Cumulative distribution functions for the completion 
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networks, all nodes that have stabilized in a state are ei- 
ther part of a stable motif, or were stabilized by the influ- 
ence of one of these motifs or by an input signal. If these 
were the only ways the state of a node can stabilize, then 
the remaining nodes would have to oscillate in the at- 
tractors. It is possible, however, (though not very likely) 
that the state of a node stabilizes without the influence 
of a stable motif. Speciflcally, this can happen when a 
node is downstream of an oscillating SCC that does not 
visit all possible states of their sub-state-space in an at- 
tractor. We discussed this kind of oscillatory behavior in 
section [II G| and called it incomplete oscillations. 

The simplest way to obtain this kind of incompletely 
oscillatory attractors is to have two nodes with self loops 
that also form a feedback loop between themselves. For 
Kauffman K — 2 networks one can show that this node 
conflguration with the choice of Boolean rules that show 
these incomplete oscillations appears in a network with 
probability 1/32A^(A^ - 1). If one adds to this that the 
downstream nodes of these node configurations need to 
have very specific update rules in order to stabilize as a 
consequence of these oscillations, it is then not surprising 
that we found only very few of these cases in our network 
ensembles. More work in this direction is needed to gen- 
eralize the reduction method to also identify these sets 
of stable nodes. 

Based on our discussion in section [lIGI one may won- 
der if unstable oscillations are also a problem for our re- 
duction method, and exactly how much of a problem they 
are. We have found that, for the ensemble of Kauffman 
critical networks, unstable oscillations are even more rare 
than incomplete oscillations. For example, the simplest 
and most probable way to have unstable oscillations (the 
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one shown in the example in Figure^ appears only with 
probability 1/128A^(A^ - 1) in a X = 2 Kauffman net- 
work. To make sure that unstable oscillations are rare, 
we searched for the SCCs that could display this behav- 
ior. For most networks (~ 90%) we have found no such 
SCCs, and even when they exist, there are only a few of 
them, as illustrated on Figure [TOl Not only that, but the 
few that have the possibility to display unstable oscilla- 
tions do not actually seem to do so (satisfying property 
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present, the asynchronous update order or of the initial state. 



(1) and (2) but not property (3) of the conditions stated 
in section |IIG| is only a sufficient condition to display 
unstable oscillations), as none of them were found by 
the attractor sampling methods in these networks. More 
work in this direction is needed to find more stringent 
conditions that identify the unstable oscillations. 

Using our network reduction method on the T-LGL 
leukemia network we were able to find the recurring sta- 
ble motifs (Figure [6]). A natural question to ask is if 
the stable motifs found during network reduction have 
any special biological significance. Interestingly, all of 
these motifs actually play a significant role in the biol- 
ogy of T-LGL leukemia: PDGFR-SlP-SPHKl-Ceramide 
(ceramide/sphingomyelin pathway) has been shown to be 
essential for T-LGL cell survival [53]; moreover, TBET 
(T-box transcription factor) and IFNG-P2 are related to 
the control of two of the main cytokines produced by 
CTLs (interleukin 2 and interferon gamma, respectively), 
whose low production is one of the characteristics of T- 
LGL leukemia [III |Ml ES] . 

An interesting observation is that these three motifs 
are directly regulated by five of the input signals of the 
T-LGL network (or almost directly in the case of TBET, 
see Figure pi) , which suggest their importance in cell-fate 
determination for GTLs. This appears to be especially 
true for the PDGFR-SlP-SPHKl-Ceramide motif, whose 
components need to always be in a specific set of states 
(PDGFR=S1P=SPHK1=0N, Ceramide=OFF) for the 
T-LGL leukemia cell fate to be possible (see Figure 11 ). 



This is consistent with and actually explains the previ- 
ous finding by Zhang et al. that an intermittent signal 
of PDGF (coupled with the sustained presence of IL15) 
was enough for T-LGL leukemia to be possible: this hap- 
pens because the intermittent signal is enough to stabilize 
the PDGFR-SlP-SPHKl-Ccramide motif in the required 
state for the T-LGL leukemia cell fate to become possi- 
ble. 

To summarize, our study showed that the novel net- 
work reduction method we propose allows us to over- 



come the limitations related to the vast state space of 
large networks by taking advantage of the stable com- 
ponents naturally present in biological networks. This is 
accomplished by transferring the complexity of the prob- 
lem from the size of the state space to the complexity of 
the network, namely, the number of cycles it has. For 
most cases, we also found that our method goes beyond 
reducing the network size and can actually predict the 
dynamical repertoire in the attractors of the system. For 
the case of the T-LGL leukemia network we found that 
the stable components identified by our method play an 
important role in the biology of T-LGL leukemia and ap- 
pear to be used as a cell-fate determination mechanism 
for CTLs. Overall, our method adds a powerful tech- 
nique to the set of tools available to infer the dynamical 
behavior of a network based on the topology and the na- 
ture of the interactions, that is flexible enough that it 
can be applied to a large variety of regulatory and signal 
transduction networks. 
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APPENDIX 

Appendix A - Proof of the conservation of attractors 
by the expanded network/network reduction method 

In the following we use V — {vi,V2, . ■ ■ ,vn) to rep- 
resent the N nodes of the Boolean network, S — 
(cri,CT2, . . . jCTat) to represent the states of these nodes, 
and F = (/i, /2, • . . , /iv) to represent the Boolean func- 
tions associated to each of these nodes. 
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We assume, for convenience, that the Boolean func- 
tions F satisfy these three properties: 

1. The /i's do not take constant values (i.e. /i 7^ 
and /, ^ 1). 

2. If fi depends on the state of node Vj, aj, then there 
must be at least a pair of states S^^^ and E^^^ with 
a J 7^ cTj , and af.' = aj,' ior k j^ j, such that 

/j(E(^^) 7^ /i(E^^^). This is equivalent to requiring 
that the Boolean derivative of fi with respect to aj 
is nonzero for at least a pair of network states [55] ■ 

3. If for any configuration of the states of a subset of 
the inputs of f^, {ai^ , cr^^ , • . ■ , ^i^. }, one has fi = 1, 
then the disjunctive normal form of fi must have 
at least one of its OR separated terms be equal to 1 
when evaluated at the state of this subset of nodes. 

The first property makes sure we have no source nodes. 
For our purposes this can be assumed without loss of 
generality, because even if that is not the case, we can 
use the reduction method of Saadatpour et al. [56] and 
remove all source nodes while preserving all attractors. 
The second property can also be assumed without any 
loss of generality; it is just a way of stating that we con- 
sider fi to depend on aj only if it explicitly depends on aj 
for at least a pair of network states. The third property 
is also general, since one can construct the respective dis- 
junctive normal form from the truth tabic of the Boolean 
function. 

Our first proposition states that the stable motifs 
found from the expanded network are such that the cor- 
responding states of these motifs are partial fixed points 
of the Boolean rules of the nodes involved. 

Proposition 1. Let M = {Vmi , Vm^ , ■ ■ ■ , Vml ) ^e a sta- 
ble motif in the expanded network representation (where 
VMi can either he a normal node, a complementary node, 
or a composite node). We denote Mstate — (cmi = bmi, 
o-ma =bm2,---,'^nii = bmi) , With 6^^. € {0, 1} OS the Cor- 
responding state of Ai in the network state E; hm =1 */ 
it is a normal node, and bm — if it is a complementary 
node. Then, for any normal node Vm or complemen- 
tary node Vjn in Ai and for any network state Tim in 
which amk — bm^yruk G {7711,7712, . . . ,7n;}, we will have 

fmj{TM) — bmj. 

Let us sketch the proof for this proposition. First, 
because a stable motif only contains a node or its com- 
plementary node, one can do a change of variables in 
the original network so that the state of the nodes 
in A^ is 1 in AA state- This simplifies the problem, 
since we now just need to show that fm{TM) = 1- 
Now, the Boolean function of node Vm ■ , fm ■ ; has the 
form /,„^ = S*! OR 52 OR • • • OR Sn, where S, = 
si AND S2 AND • • • AND si and where the s^'s are ei- 
ther a node state or its negation. Since every node Vmj 
in Ai has at least an input from other node in A^, then 
this means that one of the S^'s of fm corresponds to this 



node input. If we call Sj to the corresponding Si, then 
it must be such that all the Sk's of this Sj must be the 
state of a node in M. As a consequence Sj{T,m) = 1, 
since in the state Sm all the states of nodes in A^ is 1. 

The reverse of this proposition is also true, that is, if 
for a given set of node states updating any of the states in 
the set gives back the same state, regardless of the state 
of any node outside of the set, then this set of states 
will correspond to a set of stable motifs in the expanded 
network representation: 



Proposition 2. Let A^^ 



("■mi 



^nii 1 ^1712 ^1712 5 



be the state of a set of nodes such that 



if TiM is any network state in which cr„ 



bniyruk € 



{mi, 77i2, . . . , 777;}, then /„ . (Em) = &m, • Then (i) there 
is a set of stable motifs {Aln} in the expanded network 
representation such that each of the Ain contain only 
normal nodes or complementary nodes of the nodes whose 
state is specified in Ai state (normal nodes if 6„i^ — 1, 
and complementary nodes if &,„^_ = 0^ and in which all 
other nodes in the Aim's (if any) will be composite nodes 
made up of the normal nodes or complementary nodes 
in each Ain, and (ii) the nodes whose state is specified 
in A4 state but that are not included in the set of stable 
motifs {A^ri} will be downstream of the nodes in at least 
one of the stable motifs. 

Part of the proof for this proposition is very similar to 
the one of Proposition[l] First, one does the same change 
of variables and writes down the Boolean function of an 
arbitrary element of the nodes whose state is specified in 
Ai state- Then, from the form of Boolean function and 
since /m (Ea/) = 1, one of the parts of this Boolean 
function will be composed only of normal nodes of the 
nodes whose state is specified in A4 state, and composite 
nodes composed of these normal nodes. Finally, if one 
creates the network composed only of these normal nodes 
and composite nodes, and separates them into SCCs, one 
will find a set of source SCCs. Since sources SCCs contain 
all of its inputs (by definition) , and the elements of these 
SCCs contain only normal nodes and composite nodes 
composed of these normal nodes (by construction) , then 
these source SCCs will be the stable motifs we are looking 
for. 

For the next propositions we need to remember certain 
properties of the attractors of the asynchronous updat- 
ing scheme. For any attractor we can divide the N nodes 
into two classes: those that take the same value in all net- 
work states of ^ (i.e, either and 1), and those that take 
more than one value in the different network states of A 
(i.e, both and 1). We refer to the former as stabilized 
nodes, and to the latter as oscillating nodes. The fol- 
lowing propositions state that stabilized nodes can have 
inputs from stabilized nodes or oscillating nodes (Propo- 
sition pi), while oscillating nodes must have at least one 
oscillating node as an input (Proposition HI) . 

Proposition 3. Let A be an attractor of the Boolean 
network (V, E, F), and let S and O be the set of the sta- 
bilized and oscillating nodes in the attractor, respectively. 
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If Vs € 5, and bg is the stabilized state of node Vs, then 
one of the following two cases holds: (i) one of the OR 
separated parts of fs (ifbg = 1) or f ^ (ifbg = Oj depends 
only on the state of nodes in S , specifically, on the state 
of the nodes of S in A. If (i) is not true, then (ii) at least 
one of the OR separated parts of both fs and f ^ depend 
at least on the state of one node in O and, if it depends 
on any more states, they have to be the state of the nodes 
of S in A. 

Proposition 4. Let A be an attractor of the Boolean 
network (V, E, F), and let S and O be the set of the sta- 
bilized and oscillating nodes, respectively. If Vo ^ O then 
(i) fo O'^d fg cannot depend only on the state of nodes in 
S, (ii) neither fo nor fo can have any of their OR sepa- 
rated parts depend only on the specific state of the nodes 
of S in A (i.e, on Os ifbg = 1, or as ifbg = 0), and (Hi) 
at least one of the OR separated parts of both fo and f ^ 
depend at least on the state of one node in O and, if it 
depends on any more states, they have to be the state of 
the nodes of S in A. 

We now proceed to prove the three lemmas that will 
allow us to show that the reduction method conserves all 
attractors. In Lemma[T]we construct the set of nodes, for 
an arbitrary attractor, whose state will be identified by 
our reduction method, Sred C S. We also show that there 
is at least one stable motif composed of the corresponding 
states in the attractor of the nodes of Sred- In Lemma [2] 
we show that the network reduction of these stable motifs 
can only stabilize nodes in Sred- 

Lemma 1. Let A be an attractor of the Boolean network 
(y,E,i^), and let S and O be the set of the stabilized 
and oscillating nodes of A, respectively. There exists a 
set of nodes Sred C S such that in the expanded network 
representation of {V, ^,F) there will be at least one stable 
motif composed only of the corresponding states in A of 
the nodes of Sred, or composite nodes composed of such 
nodes. 

Proof. Without loss of generality we can do a change of 
variables so that <7s = 1 ii Vg G S. By Proposition [31 we 
can divide S into the nodes that have an OR separated 
part in their rule that depends only on the specific state 
of nodes of 5 in ^ or the nodes that have at least an OR 
separated part in their rule that depends on the state of 
a node in O. We will refer to the former as Sq and to the 
latter as Sosc- 

Let Si C So be the nodes that have at least one OR 
separated part in their rules that depends only on the 
specifc state of the nodes of Sq in A (i.e, on as, because 
of the change of variables). Let ^2 C Si be the nodes 
that have at least one OR separated part in their rules 
that depends only on the specific state of nodes of Si 
in A (note they could depend on the states of nodes in 
Sq—Si). We do this iteratively until Si^^^+i — Si^^^ and 
denote Sred — Si^^^ C Sq. Since Sred was constructed 
by first removing the nodes that required nodes in Sqsc 
to stabilize, and then removing the ones that depended 



on the previously reduced nodes, and so on, then Sred 
corresponds to the set of nodes in Sq that do not depend 
in any way on nodes of Sosc to stabilize in their state on 
A. 

We can now show that the expanded network repre- 
sentation of the states of the nodes of Sred in A has at 
least a stable motif composed only of the corresponding 
states of S in A or composite nodes composed of such 
nodes. Note that because of our change of variables, we 
only need to consider normal nodes in S and not comple- 
mentary nodes of the nodes in S. We first note that the 
construction of Sred makes sure that its nodes have at 
least one OR separated part in their rules that depends 
only on the state of nodes of Sred in A. As a consequence, 
and since the stabilized state was taken to be 1, the ex- 
panded network representation of the state of the nodes 
of Sred in A will have an input from either the nodes of 
the corresponding states of Sred in A, and/or composite 
nodes composed only of said nodes. 

From the expanded network representation we take the 
nodes of the corresponding states of Sred in A (normal 
nodes) and the composite nodes composed only of said 
nodes, and construct a network V with these normal 
nodes, composite nodes, and the edges between them. 
From the discussion in the previous paragraph, each of 
the nodes in V must have, at least, an input node. This 
means that if we divide V into SCCs there will be at least 
one sec in which the nodes have no inputs outside the 
sec itself (source SCCs). These source SCCs are stable 
motifs since the composite nodes in each of these SCCs 
have all their inputs included (the nodes in these SCCs 
have no inputs outside the SCC itself) and by construc- 
tion these SCCs contain no complementary nodes. D 

Lemma 2. Let Sred C S be the constructed set of nodes 
in Lemma uj Then (i) Sred is such that the network 
reduction of any stable motif composed only of the cor- 
responding states of Sred in A (or composite nodes com- 
posed of such nodes) can only stabilize nodes in Sred, 
and (ii) if any of the states of the nodes in Sred stabi- 
lizes, then it has to be on their corresponding state in 
A; if they do not .stabilize, then either their rule (if their 
stabilized .state in A is 1) or the negation of their rule 
(if their stabilized state is 0) will have an OR separated 
part that only depends on the .specific state of the nodes 
of Sred in A (i.e, on as if bg — 1, or as if bs = 0) that 
did not .stabilize during network reduction. 

Proof. Without loss of generality we will again do a 
change of variables so that cTs = 1 if Wj e 5. Let us 
start with the nodes in Sq — Sred (the stabilized nodes 
not in Sred), as constructed in Lemma IT] These nodes 
cannot have an OR separated part in their rules that de- 
pend only on the specific state of the nodes of Sred in A, 
since that would make them part of Sred- This means 
that setting the states of the nodes of Sred in A in their 
rules will not set the value of the rule to 1. It can also 
not make it be equal to 0, since each of them has an OR 
separated part in their rule that can be equal to 1 if the 
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state of the rest nodes they have as an input (which have 
not been set to any value) take any of the states in A. 
Hence, the nodes Sq — S^ed will not be stabilized by the 
reduction of any of the stable motifs being considered. 

Let us now look at what happens to the nodes in Sosc 
during the reduction of the stable motifs of interest. In 
Proposition [4] we showed that for every node Vq G Sqsc 
neither /o nor /^ could have one of their OR separated 
parts depend only on the specific states of the nodes of 
S in A. As a consequence, none of the nodes in Sosc will 
be stabilized by the reduction of stable motifs made up 
only of nodes in Sred C S. Since we have now shown 
than neither the nodes in 5o — S^-ed nor the nodes in Sosc 
can be stabilized by the reduction of the nodes in Sred 
in their state in A, then point (i) of the lemma has been 
proved. 

To show point (ii) , let us consider the iterative process 
involved in network reduction. First, one sets the states 
of the nodes in the chosen source SCC in the rules of all 
nodes in Sred- As a consequence, the nodes in S^ed that 
have any OR separated part of their rule depending only 
on the specific state of the nodes in the source SCC will 
stabilize in the 1 state. The rest of the nodes in S^ed 
cannot stabilize on either 1 (or they would be part of the 
previous nodes) nor on (since they still have an OR 
separated part in their rule that depends on specific the 
state of the nodes of Sred in -4 not in the source SCC). 
If one now evaluates the states of the nodes that just 
stabilized on 1 in the rest of the nodes in Sred^ and follows 
the same arguments, the result will be a new set of nodes 
that just stabilized on 1, and a set of nodes that have 
an OR separated part of their rule only depending only 
on the specific the state of the nodes of Sred in A that 
have not stabilized. Doing this iteratively until no more 
nodes stabilize one finds the desired result, that is, that 
the nodes in Sred either stabilize at their state in A or if 
they do not, then they still have one OR separated part 
that now depends only on the specific state of the nodes 
of Sred in A that this reduction did not stabilize. D 

Lemma 3. Let A he an attractor of the Boolean network 
(y, S, F), and let S and O be the set of the stabilized and 
oscillating nodes, respectively. Let Sred G S be the con- 
structed set of nodes in Lemma [7] and assume that Sred 
is empty and that O is a non empty set. Then the ex- 
panded network representation of (V, S, F) must he such 
that the normal nodes and complementary nodes of the 
elements in O, and the nodes corresponding to the state 
of the nodes of S in A must both he downstream of an 
oscillating motif that contains at least one of the nodes 
inO. 

Proof. We first show that the expanded network has, at 
least, one source SCC which has either a normal node 
or complementary node of a node in O. First, since the 
Boolean network has no source nodes, then the expanded 
network representation will also have no source nodes. 
We can then divide the expanded network into SCCs, 
with at least one of these being source SCCs. Because of 



the absence of source nodes, all nodes in the expanded 
network will be downstream of one of these source SCCs. 
Let us assume that one of these source SCCs has neither 
a normal node nor a complementary node of the elements 
in O. This means that all the inputs of these nodes are 
part of this source SCC. But, because of the construction 
of Sred and because it is taken to be empty, all nodes in S 
are downstream of a normal node or complementary node 
of the elements in O. Therefore, this source SCC must 
contain a normal or complementary node of the elements 
in O, which is a contradiction. Therefore, all source SCCs 
contain either a normal node or a complementary node 
of the elements in C As a consequence, all nodes in the 
expanded network are downstream of these source SCCs. 

Before proceeding, we first need to note a cer- 
tain property of the Boolean function with the prop- 
erties we specified at the beginning of the section. 
Without loss of generality, / can be written in the 
form / = 5i ORS-z OR ••• 0R5„, where S^ = 
si AND S2 AND • • • AND s/, and where each Sj is ei- 
ther a node state or the negation of a node state. Now, 
lets say Sk is one of the Sj's of /, then / will have one of 
its corresponding s^'s be the negation of Sk- This can be 
proved by using property 3 of the Boolean functions. 

Let Vo be any of these source SCCs that contains a 
normal or complementary node of the elements in O. 
Because of the property discussed in the previous para- 
graph, the complement of the normal nodes or comple- 
mentary nodes of Vq will also form a source SCC. These 
two SCCs have to be connected to each other, or these 
nodes would not be able to oscillate in the attractor A. 
Since they are source SCCs, this means the are actu- 
ally part of the same source SCC. This shows that Vq 
contains both the nodes and complementary nodes of ev- 
ery element it contains. Finally, since these source SCCs 
were constructed by separating the whole expanded net- 
work into SCCs, then they are the largest SCC (i.e., the 
usual definition of SCC). Hence, we have shown that 
these source SCCs are oscillating motifs. D 

The following theorem is the main result of this sec- 
tion, and it combines the results of Lemma[TJ[2] and[3j It 
shows that for every attractor in the network our reduc- 
tion method will find a corresponding quasi-attractor in 
which the state of the nodes in Sred is the same as in the 
attractor, and in which the rest of the nodes will either 
be part of an oscillating motif or downstream of it. 

Theorem 1. Let A he an attractor of the Boolean net- 
work (y, S, -F) , and let S and O he the set of the stabilized 
and oscillating nodes, respectively. Let Sred G S be the 
set of nodes constructed in Lemma [71 Then, there exists 
a set of stable motifs such that, by applying network re- 
duction, all the nodes in Sred will stabilize in their steady 
state in A, while the rest of the nodes in V will be part 
of the final reduced network. This resulting final reduced 
network will be such that, in its expanded network repre- 
sentation, all the nodes will either he part of an oscillat- 
ing motif containing at least one of the nodes in O, or be 
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downstream of an oscillating motif. 

Proof. Using Leninia[2] the network obtained after reduc- 
ing any stable motif composed only of the corresponding 
states of Sred in •A. will have a new Sred containing only 
the nodes in the previous Sred that did not stabilize. If 
one performs network reduction using the stable motif 
that necessarily exists in the new network (because of 
Lemma Rj), and does this iteratively, one will obtain a 
network where Sred is empty and where only the states 
of the nodes in the original Sred stabilized during reduc- 
tion in their state in A. By the results in Lemma [31 this 
resulting network has a set of oscillating motifs with at 
least one of the nodes in O, and with the rest of the nodes 
downstream of these oscillating motifs. D 



Appendix B - The full net^vork reduction algorithm 

In the following we describe the full network reduction 
algorithm. During the description of the algorithm we 
refer the reader to the subsections in section |TT] where the 
each of these steps are described in more detail. 



1 



For every combination of the states of the source 
nodes (nodes with no upstream components) ap- 
ply the two steps of network reduction method de- 



scribed in section II F recursively until neither of 
them can be applied anymore. 

2. Create the expanded network representation for 
each of the resulting networks (section II D ) . 



3. Search the expanded network for stable motifs (sec- 
tion II E I and oscillating components (section II G I 



4. For every separate stable motif create a copy of the 
current network. On each of the networks created 
use the states of the corresponding stable motif as 
inputs and apply the two steps of the network re- 
duction described in section |IID| recursively until 
neither of them can be applied anymore. 

5. For every oscillating component of more than two 
nodes (i.e., every oscillating component that could 
display incomplete oscillations) create a copy of the 
current network. On each of the networks created, 
the nodes in the corresponding oscillating compo- 
nent and the nodes downstream of this component 
will be marked. The marked nodes cannot be re- 
duced at any later step of the algorithm (i.e, they 
will have their state undetermined in the quasi- 
attractors that are derived from these networks). 

6. For every oscillating component of two nodes (i.e, 
only one normal node and its corresponding com- 
plementary node), check if any node downstream 
of the nodes in the oscillating motifs has a stable 
motif with no composite nodes. If none of them do, 
create one copy of the current network and mark 



the nodes in the oscillating motifs found in this step 
and the nodes downstream of them. The marked 
nodes cannot be reduced at any later step of the 
algorithm (i.e, they will have their state undeter- 
mined in the quasi-attractors that are derived from 
these networks). 

7. Repeat 2, 3, 4, 5 and 6 for each of the networks iter- 
atively until no more stable motifs are found. The 
result, a set of fixed state nodes and their stabilized 
states, and a set of nodes with undetermined states 



are the quasi-attractors (section II F I . 



8. Prune the set of quasi-attractors of duplicates (two 
quasi-attractors are the same if they have the same 
set of fixed state nodes and the same state for these 
stabihzed nodes). 



Appendix C - The logical rules of the T-LGL 
leukemia survival network 



Logical rules governing the state of the T-LGL 
leukemia survival signaling network depicted in Figure 
[5] For simplicity, the nodes states are represented by 
the node names. The Boolean rules were constructed 
based on experimental results of the corresponding 
cellular elements in healthy and leukemic cytotoxic T 
cell lymphocytes, in such a way that that the model 
reproduces the result of knockout and overexpression 
experiments. The 'NOT Apoptosis' clause in each rule 
implements the fact that in the apoptosis (cell death) 
state every node except Apoptosis is OFF. This table is 
adapted from [TTl |41] . The interested reader is referred 
to [TT] for the detailed explanation of the rules. 

fcTLA4 = TCR AND NOT Apoptosis 
/^c-j? = StimuH AND NOT (CTLA4 OR Apoptosis) 
fpDGFR = (SIP OR PDGF) AND NOT Apoptosis 
fpYN = (TCR OR IL2RB) AND NOT Apoptosis 

fcytoskeletonsignaling = FYN AND NOT ApOptOSis 

/^P^ = (CD45 OR ((TCR OR IL2RB) AND NOT 

ZAP70)) AND NOT Apoptosis 

fzAPTQ ^ LCK AND NOT (FYN OR Apoptosis) 

fGRB2 = (IL2RB OR ZAP70) AND NOT Apoptosis 

fpLCGi = (GRB2 OR PDGFR) AND NOT Apoptosis 

fRAS = (GRB2 OR PLCGl) AND NOT (GAP OR 

Apoptosis) 

fcAP ^ (RAS OR (PDGFR AND GAP)) AND NOT 

(IL15 OR IL2 OR Apoptosis) 

fhiEK ^ RAS AND NOT Apoptosis 

fERK = (MEK AND PI3K) AND NOT Apoptosis 

/p^3^ ^ (PDGFR OR RAS) AND NOT Apoptosis 

fNFKB ^ ((TPL2 OR PI3K) OR (FLIP AND TRADD 

AND lAP)) AND NOT Apoptosis 

fNFAT = PI3K AND NOT Apoptosis 

fpANTES = NFKB AND NOT Apoptosis 

/^^2 = (NFKB OR STAT3 OR NFAT) AND NOT 

(TBET OR Apoptosis) 
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fiL2RBT = (ERK AND TBET) AND NOT Apoptosis 
fiL2RB = (IL2RBT AND (IL2 OR IL15)) AND NOT 
Apoptosis 

fiL2RAT = (IL2 AND (STATS OR NFKB)) AND NOT 
Apoptosis 

fiL2RA = (IL2 AND IL2RAT) AND NOT (IL2RA OR 
Apoptosis) 

/j^^ = (IL2RA OR IL2RB OR RANTES OR IFNG) 
AND NOT (SOCS OR CD45 OR Apoptosis) 
fsocs = JAK AND NOT (IL2 OR IL15 OR Apoptosis) 
fsTATs = JAK AND NOT Apoptosis 
/P27 = STAT3 AND NOT Apoptosis 
fprouferauon = STAT3 AND NOT (P27 OR Apoptosis) 
Itbet = (JAK OR TBET) AND NOT Apoptosis 
fcREB = (ERK AND IFNG) AND NOT Apoptosis 
fiFNGT = (TBET OR STAT3 OR NEAT) AND NOT 
Apoptosis 

fiFNG = ((IL2 OR IL15 OR Stimuli) AND IFNGT) 
AND NOT (SMAD OR P2 OR Apoptosis) 
/P2 = (IFNG OR P2) AND NOT (Stiniuli2 OR Apop- 
tosis) 

fazMB = ((CREB AND IFNG) OR TBET) AND NOT 
Apoptosis 

/tpl2 = (TAX OR (PI3K AND TNF)) AND NOT 
Apoptosis 

fj.j^p = NFKB AND NOT Apoptosis 
Jtradd = TNF AND NOT (lAP OR A20 OR Apopto- 
sis) 

fpasL = (STAT3 OR NFKB OR NFAT OR ERK) AND 
NOT Apoptosis 

/p^,-r = NFKB AND NOT Apoptosis 
fp^^ = (FasT AND FasL) AND NOT (sFas OR Apop- 
tosis) 

fsFas = FasT AND SIP AND NOT Apoptosis 
fceram^de = Fas AND NOT (SIP OR Apoptosis) 



/disc = (FasT AND ((Fas AND IL2) OR Ceramide 

OR (Fas AND NOT FLIP))) AND NOT Apoptosis 

fcaspase = ((((TRADD OR GZMB) AND BID) AND 

NOT lAP) OR DISC) AND NOT Apoptosis 

fp^jp = (NFKB OR (CREB AND IFNG)) AND NOT 

(DISC OR Apoptosis) 

/^2o = NFKB AND NOT Apoptosis 

Ibid ^ (Caspase OR GZMB) AND NOT (BclxL OR 

MCLl OR Apoptosis) 

fiAP ^ NFKB AND NOT (BID OR Apoptosis) 

JbcIxL = (NFKB OR STAT3) AND NOT (BID OR 

GZMB OR DISC OR Apoptosis) 

Jmcli = (IL2RB AND STAT3 AND NFKB AND 

PI3K) AND NOT (DISC OR Apoptosis) 

f Apoptosis = Caspase OR Apoptosis 

Igpcr = SIP AND NOT Apoptosis 

fsMAD = GPCR AND NOT Apoptosis 

fsPHKi = PDGFR AND NOT Apoptosis 

/sip = SPHKl AND NOT (Ceramide OR Apoptosis) 

Appendix D - The attractors of T-LGL leukemia 
survival network 

In Table |l] we show the state of the nodes for all possi- 
ble combinations of input signals in the presence of anti- 
gen (Stimuli=ON). We do not show the Apoptosis=ON 
attractor, in which all nodes except Apoptosis are inac- 
tive, since it is always a possibility. For simplicity, we 
only show which nodes oscillate and which of them sta- 
bilize in a steady state (i.e., the quasi-attractor) and not 
the actual attractor, which would include all the network 
states that the nodes that oscillate can visit along with 
the transitions between these states. The signal com- 
binations not shown in the table (i.e. CD45=OFF and 
IL15=0FF, with any other value for the other input sig- 
nals) have apoptosis as their only attractor. 
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TABLE I: The attractors of T-LGL leukemia survival network. This table shows the state of the nodes for all 
possible combinations of input signals in the presence of antigen (Stimuli=ON). 
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